Assigned Reading: Geocomputation with R Chapter 2.

Sacramento Mountains, New Mexico

Sacramento Mountains, New Mexico

What is a Geographic Information System?

What is a GIS

A Visualization of a GIS, by Anne Sexton

A Visualization of a GIS, by Anne Sexton

A Brief History

The Broad Street Pump, A Cholera Outbreak in London, By: John Snow and digitized by National Geographic

The Broad Street Pump, A Cholera Outbreak in London, By: John Snow and digitized by National Geographic

https://www.theguardian.com/news/datablog/2013/mar/15/john-snow-cholera-map

What is Spatial Data Science?

Why Use R as a GIS and for Spatial Data Science ?

Geodesy

Earth is not a perfect sphere

Blue Marble 2012, by Suomi NPP

Blue Marble 2012, by Suomi NPP

The Earth can be represented as an Ellipsoid

Ellipsoids

Ellipsoids

However, the Earth’s surface is not smooth - Geoid

Left: Geoid Cross Section. Right: IGCM GeoidLeft: Geoid Cross Section. Right: IGCM Geoid

Left: Geoid Cross Section. Right: IGCM Geoid

By International Centre for Global Earth Models (ICGEM) - http://icgem.gfz-potsdam.de/vis3d/longtime / Ince, E. S., Barthelmes, F., Reißland, S., Elger, K., Förste, C., Flechtner, F., Schuh, H. (2019): ICGEM – 15 years of successful collection and distribution of global gravitational models, associated services and future plans. - Earth System Science Data, 11, pp. 647-674,DOI: http://doi.org/10.5194/essd-11-647-2019., CC BY 4.0, https://commons.wikimedia.org/w/index.php?curid=81462823

Geodetic Datums

Other locations are then measured in relation to the control points.

Number one City Datum. by: Cosmo1976

Number one City Datum. by: Cosmo1976

Control Points Number one City Datum

Control Points Number one City Datum

Geodetic Coordinates

“The geodetic (or geographic) latitude is the angle between the equatorial plane and the normal (vertical) to the ellipsoid surface at the considered point.” - Proj

Angles on Ellipsoid and Geodetic Coordinates. by: Peter MercatorAngles on Ellipsoid and Geodetic Coordinates. by: Peter Mercator

Angles on Ellipsoid and Geodetic Coordinates. by: Peter Mercator

By Peter Mercator - Own work, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=17717979

Coordinate Notation

\[ \text{Decimal Degrees} = \text{Degrees } + \frac{\text{Minutes}}{60} + \frac{\text{Seconds}}{3600} \] \[ \text{Latitude of Tech} = 42°03'27.7" N \]

\[ 42 + (\frac{03'}{60}) + (\frac{27.7"}{3600}) \]

\[ 42 + 0.05 + 0.007694 = 42.05759 \text{ Decimal Degrees} \]

UTM - Common for field work, planar measurements of meters.

UTM Zones of the Conterminous USA. data by ESRI

UTM Zones of the Conterminous USA. data by ESRI

Coordinate Reference System

Geographic Coordinate Reference System. by: Anna Krystalli

Geographic Coordinate Reference System. by: Anna Krystalli

Problems with working on flat surfaces

3D to 2d does not work well. But we have worked around this. sorta

Orange Peel. by: Nathan Belz

Orange Peel. by: Nathan Belz

Geographic Coordinate System

Projected Coordinate System

Major Map Projections

Many different ways to create two dimensional maps – thousands of projections. None are perfect. There are four main types of projections, each with their own strengths and examples. One of the most common examples of each of these is included in the table below.

Map Projections. by: Daniel Strebe

Map Projections. by: Daniel Strebe

Geodesy Takeaways

Unless you really focus on GIS, hardly anything I have said this lecture will matter to you.

An Introduction to Geographic Data Models

Vector data tends to only include features of interest, e.g. bodies of water; whereas a Raster will include, explicitly, the absence of features (e.g. both water and terrestrial areas).

Vector and Raster Data Models

Vector and Raster Data Models

Vector Data in R

Simple Features - Standards

Simple Features 1 - Attributes

TAXON DBH HEIGHT
Robinia pseudoacacia 40 24
Quercus alba 32 21

Simple Features 2 - Coordinates form a Point

Simple Features 3 - Points are/form a SF Geometry (sfg’s)

Simple Features 4 - A Geometry can be combined to form Geometries (sfg’s)

Simple Features 5 - Geometries and Spatial Information form an SF Collection (sfc)

Wherein precision refers to the precision of the coordinates forming the geometries. the bbox is a rectangle which includes all pairs of coordinates in the geometry. N_empty, denotes whether the geometry is missing coordinates in a position.

Simple Features 6 - Recapped.

TAXON DBH HEIGHT LONG LAT geometry
Robinia pseudoacacia 40 24 -87.67607 42.05663 POINT (-87.67607 42.05663)
Quercus alba 32 21 -87.67399 42.05715 POINT (-87.67399 42.05715) 
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -87.67607 ymin: 42.05663 xmax: -87.67607 ymax: 42.05663
Geodetic CRS:  WGS 84
Precision:     2 
sfc_POINT of length 1; first list element:  'XY' num [1:2] -87.7 42.1

sp

sp - History

sp - Spatial Classes for Topology

sp - Structures of The Spatial* Class - Topology Only

## Formal class 'SpatialPoints' [package "sp"] with 3 slots
##   ..@ coords     : num [1:15, 1:2] 0.49 0.38 0.27 0.11 0.64 0.03 0.79 0.38 0.28 0.16 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "xc" "yc"
##   ..@ bbox       : num [1:2, 1:2] 0.03 0.09 0.93 1.78
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "xc" "yc"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

sp - Attributes - DataFrame

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  15 obs. of  2 variables:
##   .. ..$ temperature: num [1:15] 7.02 4.49 5.4 4.71 4.33 4.73 4.34 5.4 3.92 5.49 ...
##   .. ..$ aspect     : int [1:15] 16 11 5 17 26 22 38 42 8 35 ...
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:15, 1:2] 0.49 0.38 0.27 0.11 0.64 0.03 0.79 0.38 0.28 0.16 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "xc" "yc"
##   ..@ bbox       : num [1:2, 1:2] 0.03 0.09 0.93 1.78
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "xc" "yc"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA
temperature aspect
7.02 16
4.49 11
5.40 5
4.71 17
4.33 26
4.73 22

We see here that SF contains all of the spatial information in a Simple Feature Collection which can be tidily tucked away in the geometry list column, in SP the analogue to the SFC is the object itself.

SP is mostly geometry with a data slot, SF is mostly data frame with a geometry column.

sp Overview of a Spatial*DataFrame

sp - Recapped.

Raster data in R

Western Plant Predictors. Note each cell is the extent of a single raster tile, each tile contains cells.

Western Plant Predictors. Note each cell is the extent of a single raster tile, each tile contains cells.

Raster Components

Note that both cell sizes, and the resolution of the values in cells can, within reason, be converted to finer and coarser resolution. We will discuss these types of calculations next class.

Example Raster 1 Create Frame

Rasters tend to confuse people. We are going to create our own example here before we talk about them much more.

Fairy shrimp live in bodies of water which dry out in late Spring, and refill in early Fall. We seek to determine where suitable habitat for fairy shrimp are in the Great Basin. In order to do so, we will use satellite imagery to detect valleys which fill with water in Fall, and dry out in late Spring.

Suprise Valley by: John Glen

Suprise Valley by: John Glen 

empty_raster <- raster(
  
  # rasters have 4 bounding edges
  # Here we define each 'corner' of the raster'
    xmn = 697129.7, 
    xmx = 811775.7,
    ymn = 4388466,
    ymx = 4502382,
    
  # Here we set the number of cells 118*118
    nrows = 118, 
    ncols = 118,
    
    # we do NOT manually specify the cell size here.
    
    crs = "+proj=utm +zone=10 +datum=WGS84", 
    # set the rasters Coordinate Reference System
)

In the code above, we are going to define all four of the essential components of a raster. Clearly, we are defining three explicitly (Bounding Coordinates, Dimensions, CRS), and the remaining one implicitly (Cell resolution).

We can imagine that the raster we are currently creating looks like this. A frame without content. We can see where the bounding coordinates are,

Example Raster 2 Set Values

raster_matrix <- matrix(rast_vals_num, # # fill matrix with values, 
                       nrow = empty_raster@nrows, # create matrix same dimensions.
                       ncol = empty_raster@ncols, # create matrix same dimensions
#ensure filling of matrix goes from upper left to lower right.
                       byrow = T) 

raster_matrix[90:118,112:118] <- 1 # fix the clipping image at edge.
Example Matrix showing values underlaying a raster layer.
0 1 1 0 1 0 0 1
1 0 1 0 0 1 1 0
1 0 0 1 1 0 1 0
0 1 1 0 1 0 0 1
1 0 1 0 0 1 1 0
1 0 0 1 1 0 1 0

I think it is easiest to imagine that the values within a raster are in the shape of a matrix. Hence, using the code above we could take a vector of values, and bend them, so that each position within the vector matches up with the beginning of a new row.

example_raster_dec <- setValues(empty_raster, raster_matrix)

Example Raster 3

The width of each raster cell is: 971.57627 meters
The height of each raster cell is: 965.38983 meters
This example_raster_dec contains: 13924 elements

Raster Package - in R!

Attributes

Raster Layer

A single raster layer

Raster Stack

In general these are rasters which have been classified and we want to extract values from or run calculations with. Now what is great about layers, is that we can do one big thing bricks cannot do, we can load in many layers from many files and create a stack of attributes we are interested in studying on the fly.

The main use of Raster Stacks is to hold layers of similar themes.

Two usage examples:

  1. Each layer is the same variable from a different time.
  • e.g. mean monthly temperature from January -> December (12 layers per stack)
  1. Each layer is a different variable (theme) in an analysis.

  • e.g. yearly mean temperature, mean precipitation, etc.

  • Do not need be held in memory

Raster Brick

  • Multiple bands of imagery held in the same file
  • The sensors of a camera, e.g. Red, Green, and Blue
  • Used for performing image classification to produce data raster layers
  • Now bricks of over 100 sensor bands exist…

As I have mentioned rasters are often generated from satellite imagery; we will discuss rasters which are not developed this way next lecture. Historically most pictures are imaged via the use of three bands. These having spectral values of Red, Green, and Blue (RGB) associated with them. For example our .tif file, is composed of three bands. Note that a Rasterbrick is most often used for loading in these types of imaging data, which can then be processed to form a more typical ‘raster’ dataset.

- do not need to be held in memory,

Terra

Cartography

data("us_states", package = "spData")
us_states <- st_as_sf(us_states)

north_carolina <- read_sf(system.file("shape/nc.shp", package = "sf"))
ggplot(north_carolina) +
  geom_sf() + 
  theme_bw() # we will use this theme for class.

ggplot(north_carolina) +
  geom_sf(aes(fill = BIR79))  + # fill the interior of polygons by a variable
  theme_bw()

ggplot(north_carolina) +
  geom_sf(aes(color = BIR79),  # color the borders of polygons by a variable
          lwd = 1 # just making the borders thicker.
          )  +
  scale_color_viridis_c(option = "plasma", trans = "sqrt") +
  # just using a very colourful color scheme so you can see it. 
  theme_bw()

ggplot(north_carolina) +
  geom_sf(aes(color = BIR79),
          fill = NA  # set to 'NA' to 'remove' the interior of polygons
          ) +
  scale_color_viridis_c(option = "plasma", trans = "sqrt") +
  theme_bw()

ggplot(north_carolina) +
  geom_sf(aes(fill = BIR79),
          color = NA  # set to 'NA' to 'remove' the interior of polygons
          ) +
  scale_color_viridis_c(option = "plasma", trans = "sqrt") +
  theme_bw()

ggplot(north_carolina) +
  geom_sf(aes(fill = BIR79, 
              color = FIPS)
          ) + # both fill and color
  guides(color = 'none') + # removed the categorical legend (is too big!)
  theme_bw()

ggplot() + # leave empty
  geom_sf(data = north_carolina, # put data here
          aes(fill = BIR79, 
              color = FIPS)
          ) + # both fill and color
  guides(color = 'none') + # removed the categorical legend (is too big!)
  theme_bw()

ggplot(us_states) + # two data sets.
  geom_sf() +
  geom_sf(data = north_carolina, 
          fill = 'purple', 
          color = NA # hash this line out in your own time to see the difference
          ) +
  theme_bw()

ggplot(north_carolina) + # oh no we drew over North Carolina!
  geom_sf(fill = 'purple') +
  geom_sf(data = us_states) +
  theme_bw()

ggplot() + # two data sets.
  geom_sf(data = us_states, fill = NA) +
  geom_sf(data = north_carolina, fill = 'purple') +
  theme_bw()

ggplot() +
  geom_sf(data = us_states) +
  geom_sf(data = north_carolina) +
  coord_sf(crs = 4267, # we can convert each dataset to the same CRS 
           datum = NA # we can remove the grid lines/graticules from plot
           )

bound <- st_bbox(north_carolina) # retrieve the bbox from the sfc list column

ggplot() +
  geom_sf(data = us_states) +
  geom_sf(data = north_carolina, aes(fill = BIR79)) +
  coord_sf( # use the bbox to 'crop' the extent of the map
           xlim = c(bound[1], bound[3]), 
           ylim = c(bound[2], bound[4])
           ) +
  theme_bw()

Assignments for the Duration of the Spatial Data Science Module:

For next Lab:

Please install these packages (if you have not done so already):

# install.packages("sf", "raster", "terra", "sp", "tmap", "leaflet", "ggmap")
# optional packages for using parallel processing at a step
# (in our example it won't actually speed up the analyses 
# they actually may slow them down!)
# install.packages('snow','parallel')

Download the labs .R script from the course website.

If you are interested in how Drone/Lidar data are collected please check out the ‘RMBL Spatial Data Science Webinar Series’:

For next Lecture: Assigned Reading: Chapter 3 of Spatial Data Science

Future Bonus SDS Office Hours Wednesday Night at 5:00 - 6:00.

Notes on this Lab My lecture notes are in the R script I used to generate all of the novel figures for this presentation. Likewise this presentation is an .HTML file and can be launched from your computer (it was rendered directly from R using the script).

Works Cited

Krystalli, A. https://annakrystalli.me/intro-r-gis/gis.html Accessed 01.20.2022

https://geocompr.robinlovelace.net/spatial-class.html Accessed 01.09.2022

Hijman, R. 05.12.2019 ‘The raster Package’

https://rspatial.org/raster/RasterPackage.pdf Accessed 01.09.2022

https://www.neonscience.org/resources/learning-hub/tutorials/dc-multiband-rasters-r Accessed 01.19.2022

Pebesma, E. https://r-spatial.github.io/sf/articles/sf1.html Accessed 01.10.2022

https://proj.org/operations/conversions/geoc.html Accessed 01.20.2022.

https://rspatial.org/raster/spatial/8-rastermanip.html Accessed 01.11.2022

https://en.wikipedia.org/wiki/Open_Source_Geospatial_Foundation Accessed 01.09.2022

NOAA. What is geodesy? National Ocean Service website, https://oceanservice.noaa.gov/facts/geodesy.html, 1/25/17.

Geocomputation with R. Lovelace, R., Nowosad, J., Muenchow, J. 2022-01-25.

Packages Cited:

c("raster", "sp", "sf", "tidyverse", "terra") %>%
  map(citation) %>%
  print(style = "text", na.print = '')
[[1]]
Hijmans R (2022). _raster: Geographic Data Analysis and Modeling_. R
package version 3.5-15, <URL:
https://CRAN.R-project.org/package=raster>.

[[2]]
Pebesma EJ, Bivand RS (2005). "Classes and methods for spatial data in
R." _R News_, *5*(2), 9-13. <URL:
https://CRAN.R-project.org/doc/Rnews/>.

Bivand RS, Pebesma E, Gomez-Rubio V (2013). _Applied spatial data
analysis with R, Second edition_. Springer, NY. <URL:
https://asdar-book.org/>.

[[3]]
Pebesma E (2018). "Simple Features for R: Standardized Support for
Spatial Vector Data." _The R Journal_, *10*(1), 439-446. doi:
10.32614/RJ-2018-009 (URL: https://doi.org/10.32614/RJ-2018-009), <URL:
https://doi.org/10.32614/RJ-2018-009>.

[[4]]
Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R,
Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E,
Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi
K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the
tidyverse." _Journal of Open Source Software_, *4*(43), 1686. doi:
10.21105/joss.01686 (URL: https://doi.org/10.21105/joss.01686).

[[5]]
Hijmans R (2022). _terra: Spatial Data Analysis_. R package version
1.5-18, <URL: https://rspatial.org/terra/>.